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A simple analytic model predicted Hall thruster channel erosion based on thruster geometry, operating conditions, 
and magnetic field configuration. This model relied on a one-dimensional representation of the plasma with a fixed 
ionization fraction and variable ion energies based on the magnetic field distribution. Sputtering was modeled as the 
result of elastic scattering of ions by neutrals within the channel. Not all scattered ions and neutrals were assumed to 
reach the channel walls as a result of additional subsequent scattering events. Incorporating this phenomenon 
resulted in a greater predicted decrease in erosion rate with time than predicted based only on geometric effects. 
Results from this model were compared to SPT-100 experimental erosion data. 


Nomenclature 

A - area 

B - magnetic field 

E - electric field 

e - electric charge 

f( (p) - incident angle distribution function 

g( e) - incident energy distribution function 

g - relative velocity 

I - current 

L - length of computational cell 

m - mass 

n - number density 

N - number flow rate 

p - pressure 

p c - probability of primary scattering collision 

p c2 - probability of secondary deflection collision 

q - ion charge 

r - radius 

S y - sputter yield 

S y - normalized sputter yield 

U - diffusion velocity 

u - velocity 

V - potential, volume 

Z - collision rate 

Az - length of wall segment 

a - random number (-1,1) (magnitude of post-collision axial velocity) 

8 - energy of incident particle 

6 - random number (0,27i) (post-collision angle from radial axis) 

A - mean free path 

v - collision frequency 

cr - collision cross section 

(p - incident angle with a wall 
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Subscripts 

cm - center of mass 

cross - channel cross-section 

e - electron 

erode - distance eroded 

i - ion 

n - neutral atom 

new - radius before erosion 

old - radius after erosion 

r - radial 

wall - of the wall 

z - axial 

£ - incident particle energy 

6 - azimuthal 

(p - incident angle of particle 


I. Introduction 

Electric propulsion systems utilizing Hall thrusters are being considered for a diverse range of missions requiring 
specific impulses ranging from 1000 to over 4000 seconds in support of NASA’s current vision for space 
exploration. 1 ' 3 One critical aspect of Hall thruster propulsion systems required to support these applications is 
the operational lifetime of the engine. The primary life-limiting mechanism of magnetic-layer Hall thrusters is 
erosion of the ceramic discharge chamber by energetic ions contained within the thruster. The rate of discharge 
chamber erosion is substantially affected by discharge voltage, which varies from approximately 100 Volts to over 
1000 Volts over the aforementioned range of specific impulse. Experimentally evaluating thruster lifetime at each 
power level and discharge voltage of interest is extremely expensive and time consuming due to the long thruster 
lifetimes (typically on the order of 10,000 hours) and facility requirements for such tests. An analytic tool capable of 
quickly simulating the Hall thruster erosion process based on operating condition could provide an estimate of 
lifetime for a fraction of the cost. Such a tool would not negate the need for long duration qualification tests, but 
would enhance NASA’s ability to evaluate the suitability of Hall thruster technology for applications of interest. 

Previous investigators have developed analytic tools to predict Hall thruster lifetime. 4 ' 7 Most have been semi- 
empirical in nature, require experimental data as an input, and have indeterminate applicability at operating 
conditions substantially different from those used to determine their empirical relationships. In order to address 
these shortcomings, a simple physics-based analytic life predictive capability was sought based on thruster 
configuration, thruster operating conditions, and known sputter yields. During the present investigation an analytic 
model that predicts time-dependent erosion rates and channel profiles, based on thruster geometry, xenon ion sputter 
yields for boron nitride, anode mass flow rate, and the radial magnetic field distribution, was developed. This paper 
describes the simplified analytic approach, discusses the assumptions made in developing the model, compares 
calculated results with experimental data, and suggests areas for future improvement. 


II. Approach 

Hall thruster erosion results from the interaction of energetic atoms and ions within the discharge chamber with 
the walls of the discharge chamber. This interaction is a sputtering process that occurs at the surface of the ceramic 
discharge chamber as a result of the impact of heavy particles (atoms and ions) with sufficient kinetic energy to 
overcome the inter-atomic coupling forces binding together the constituents of the ceramic channel walls. The rate 
of sputtering is determined by the rate at which heavy particles impact the wall, the angles of incidence, and the 
heavy particles’ energies. Therefore, in order to calculate the time-dependent Hall thruster discharge chamber 
erosion rate it is necessary to determine the collision rate with the walls, the angle of incidence for particle-wall 
collisions, and the energy of these collisions. A simple physics based approach was adopted for calculation of each. 
This approach relies on a one-dimensional model for the plasma parameters. Erosion was modeled as the result of 
elastically scattered ions and atoms impacting the channel at the same axial location they were scattered. The angle 
of incidence for collisions with the wall was calculated assuming the scattering process was isotropic. Inputs to the 
model include channel geometry, cathode location, discharge voltage, magnetic field strength distribution, and mass 
flow rate. 
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This approach was initially applied to the analysis of erosion of a Hall thruster with the discharge voltage and 
thruster geometry shown in table 1. The cathode was assumed to be located 0.024 m downstream of the channel exit 
plane. These values were chosen to permit comparison with existing erosion data. 8 

TABLE l.—HALL THRUSTER MODEL INPUTS 


Outer diameter 

0.100 m 

Inner diameter 

0.069 m 

Channel length 

0.024 m 

Discharge voltage 

300 V 


The determination of ion energies, based on position, discharge voltage, and magnetic field strength distribution 
began with the electron momentum equation. 9 


d_ 

dt 


(. m e n e u e ) + V • m e n e {u e u e 


-U e U e ) = -en e (E+u e xB)~ Vp e 


-m e n e vJU e -U i) - m e n e v en {U e -U n ) 


( 1 ) 


where m is mass, n is number density, u is drift velocity, U is diffusion velocity, e is the electric charge, E is the 
electric field, B is the magnetic field, p is pressure, and v is the collision frequency. The subscripts e , i, and n 
represent electrons, ions, and neutrals respectively. This equation was simplified considerably. Steady flow was 
assumed eliminating the time derivative. The diffusion velocities were assumed small compared to the drift 
velocities. The electron drift velocity square term and pressure gradient were both assumed to be negligibly small. 
Based on these assumptions the momentum equation was reduced to 

E = —u e x B (2) 


In the present analysis, the electron drift velocity was assumed to be independent of the magnetic field, which 
results in the electric field being proportional to the magnetic field. Based on this relationship between the electric 
and the magnetic fields, the potential was obtained through an integration of the electric field normalized by the 
discharge voltage. The potential profile determined in this fashion for the SPT-100 is shown in figure 1. 



Axial distance from thruster exit, mm 


Figure 1 . — Potential distribution versus axial position 

Based on this potential distribution ion velocities were then determined by integrating the force described by 
Ohm's law. 

^rn t uf = J q(E z + u e xB r )dz (3) 
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where q is the charge of the ion, which was assumed to be all singly charged. Since the azimuthal current of the ions 
was assumed to be negligible, the ion velocities were found from the potential drop, AV. 


u ; = 


2q 


m, 


AV 


(4) 


Neutral and ion densities were estimated based on the channel cross-section, the calculated ion velocity, total mass 
flow rate, and an assumption of 90% propellant utilization. While this assumption is not suitable upstream of the 
ionization region, the channel erosion process occurs near the thruster exit subsequent to ion creation. The resulting 
ion velocity and number densities are shown in figure 2. A constant neutral density of 4.9x1 0 18 nT 3 was used 
throughout the computational domain based on conservation of mass. These values are based on an assumed anode 
mass flow rate of 4.9 mg/s. 
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Figure 2 . — Ion density and velocity versus axial position 


The flux of atoms and ions to the walls of the discharge chambers was assumed to be the result of elastic 
scattering. This process has been shown to be important at the lower densities found in the plume, 10 and permitted a 
simplified one-dimensional approach. The probability of an ion colliding with a neutral atom was calculated based 
on the following collision cross section, 11 

a in =1.69xl0~ 18 -2.77xl0~ 19 log 10 (g) m 2 (5) 

The collision rate was then calculated: 

Z in = n i n n°in8 ( 6 ) 

where the relative velocity, g , is approximated by the ion velocity. The cross section was doubled in order to 
account for momentum exchange and charge exchange collisions, both of which are important to consider in finding 
the number flux of particles to the walls. The collision rate per ion was determined by dividing the total collision 
rate by the ion number density. This rate was then multiplied by the average residence time of an ion within a 
computational cell, namely At = L/u t where L is the length of the cell. The resulting probability, and thus the 
average fraction of ions which undergo a collision in a cell, was then determined: 

p c =n n [3.38x!0~ 18 -5.53xl0~ 19 log 10 (u^L (7) 
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A plot of the collision probability versus axial distance from the anode is shown in figure 3. The probability after 
4000 hours, affected by changes in the geometry and therefore number densities, is also shown. 



Figure 3. — Fraction of ions scattered versus axial position 


The directionality of the scattered ions was determined by assuming an isotropic process. The center- of-mass 
velocities as well as the magnitude of the relative velocity were both conserved. However, the relative velocity 
vector was randomly oriented after collision. To achieve this, a random fraction of the axial component was 
retained. The components of the post-collision velocity vector for the ions, taking the center-of-mass velocity into 
account, were 12 


* Win I I 

U z= U cm„ + 

m i + m n 

-|g| • a / 1 — CL • COS 6 


r m l + m n 


where the center-of-mass velocities are given by 


( 8 ) 

( 9 ) 


m, m„ 

U cm= f U i + U n 

Wi + Wn m i+ m n 


( 10 ) 


All pre-collision velocities were assumed to be axial, thus u cmr = 0. For each ion collision, the corresponding neutral 
atom was also tracked. The neutrals have post-collision velocity components that conserve momentum and the 
center-of-mass velocity of the two-particle system. 

From the post-collision velocity components, the resulting scattering angles were obtained. The tangent of the 
scattering angle was obtained by dividing the radial component of the post-collision velocity by the axial 
component. The compliment of this angle is the angle of incidence with the walls, <p, therefore: 


tan (p = 


m n 


+ m, 


M- a 


m„ 


niz + m, 


|g| • a/ 1-oc cosO 


( 11 ) 
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The distributions of the scattering angle for various positions along the channel are shown in figure 4. These 
distributions were obtained by sampling over a large number of trials and normalizing. 



The current to the discharge chamber walls was determined from the fraction of ions that undergo scattering 
collisions. Current to the outer discharge chamber wall was the result of ions with a positive post-collision radial 
velocity component. The inner wall current was the result of ions with a negative post-collision velocity component. 
Both the scattered ions and the corresponding scattered neutrals were assumed to hit the walls in the same axial 
location they were scattered. Conservation of momentum insured that for each ion that hit the outer wall, the 
corresponding neutral impacted the inner wall at the same axial location. The current to each segment of wall was 
determined based on the plasma current at that axial location. The plasma current was found by multiplying the 
plasma density by the cross-sectional area. The current to the walls is then found by considering the fraction of ions 
experiencing scattering collisions at that axial location, i.e., 

I wall ~ i) A c W ss Pc ( 12 ) 


where A cross is the cross-sectional area of the channel. 

An additional factor was taken into account when calculating the current to the walls. As the walls erode, the 
radial distance particles travel before impacting the wall increased. It was assumed there was an additional chance 
that scattered particles experienced another collision and were scattered away from the walls before impacting the 
wall. This was approximated by assuming a mean free path (in this case, X = 4 mm) and comparing it to the extra 
radial distance caused by erosion, r erode . All of these additional scattering collisions were assumed to deflect the 
particles away from the walls for simplicity. The probability that a particle was deflected away was then 


Pel 



^ erode 

X _ 


(13) 


The calculated rate of ion collisions for both the inner and outer discharge chamber walls are shown in figure 5 at 
the beginning of life and after 4000 hours. The calculated number flow rate of neutrals to the walls is similar to that 
of the ions. 
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Figure 5. — Number flow rate to the walls 


The energy and incident angle for the ion and neutral flows shown above were then determined. The incident 
angle distribution obtained by binning a large number of random samples, is shown in figure 6. Only the 
distributions of a few wall segments are shown, though the others follow the general trend seen in the figure. The 
distribution of neutral atom incident angles was identical. 



Figure 6. — The distribution of incident angle of ions impacting the outer wall 


The energy distribution used is shown in figure 7. For the ions, the pre-scattered energies were used. This 
approximated the acceleration the ion experiences after scattering from a further potential drop. Since the pre- 
scattered energy was used, the energies were discretized as per the domain, and appeared as delta functions on the 
distribution plot. For the neutrals, the post-collision velocities were used to determine their impact energy; therefore 
they had a distribution of energies. Both figures 6 and 7 show the result along the outer wall. The inner wall 
distributions are similar. 
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Figure 7. — The distribution of particle energy impacting the outer wall 


Figures 8 and 9 show the data of boron nitride sputter yields as a function of incident angle and energy. 13 
Corresponding curve fits are also shown. Figure 8 displays the angular dependence of the normalized sputter yield. 
Figure 9 shows the sputter yield versus ion energy at a normal incidence angle. 



Figure 8. — The angular dependence of the normalized sputter yield 
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Figure 9 . — The BN sputter yield versus ion impact energy ( normal incidence ) 


A threshold energy of 50 eV was assumed. The curve fit equations are given below. 

S y(p =-1.89xl0~ 7 (p 4 +2.04xl0~ 5 (p 3 -3.77xl0~ 4 (p 2 + 1.85 x 10~ 3 cp+ 0.426 (14) 

S ye = e(0.0346 Ine- 0.136) (15) 

where (p in this case was the energy of the incident ion. The sputter yield at a given wall location was found by 

multiplying the angle distribution, f( (p), by the normalized sputter yield angular dependence by the energy 
distribution, g(e ), and the sputter yield angular dependence. 

S y =f((p)S y(p -g(e)S y£ (16) 

The sputter yield multiplied by the flow to the walls was used to estimate the volumetric erosion rate at each 
location along the discharge chamber wall. 

dV 

- A ~ = S y N (17) 

dt y 

where V represents the volume of wall material eroded away. The erosion rate was multiplied by the time step to 
determine how much of the wall was eroded for that time interval. The new radii at a particular axial location was 
found geometrically 


r 2 +■ 
r old - 


v 

71 Az 


(18) 


where r new was the updated radius after the specified erosion, r old was the previous radius, and Az was the axial length 
of the wall segment under consideration. The positive sign was used for the outer wall and the negative sign was 
used for the inner wall. The resulting change in geometry was taken into account for the subsequent time step 
altering parameters such as number density and the probability of secondary scattering away from the walls. 
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III. Results 

Based on this analytic model, erosion profiles were calculated for an SPT-100 thruster geometry, an applied 
voltage of 300 Volts, and an anode mass flow rate of 4.9 mg/s. The calculated erosion profiles obtained for the inner 
and outer walls are shown in figure 10 and 11 along with measured values 8 at various time intervals ranging from 
160 hours of operation to 4000 hours of operation. The results, in general show good agreement with the 
experimental data. This good agreement was obtained by artificially increasing the neutral number density by a 
factor of two in order to match the measured beginning of life erosion rate. Increasing the neutral density by a factor 
of two does not seem unreasonable based on a radial gradient in neutral density and less than 90% propellant 
ionization in the region of the channel adjacent to the channel wall. 



Figure 10. — The erosion profiles for the inner wall 



Figure 11. — The erosion profiles for the outer wall 


The model did accurately predict that channel erosion begins approximately 9 mm upstream of the exit, where 
calculated ion energies exceed the 50eV sputter threshold. Substantial uncertainty exists in the xenon-boron nitride 
sputter yield at these low energies, potentially resulting in the under-prediction of erosion in this region of the 
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channel. This discrepancy may also be the result of the model’s failure to consider ions scattered from other axial 
locations. 

Closer to the thruster exit, the model initially over-predicts the erosion rate and under-predicts erosion for times 
longer than 2400 hours. The predicted change in erosion rate is due to geometry changes and the possibility of 
scattered ions and neutrals being scattered a second time before impact with the wall. While there is disagreement 
between the measured wall profiles and those predicted by the model, the overall volumetric erosion rate compares 
favorably with experimental data as shown in figure 12. The effect of secondary scattering was required in addition 
to the geometry effects in order to accurately predict the change in erosion rate over time. 



IV. Conclusions and Recommendations 

A simple analytic model was demonstrated that predicts Hall thruster channel erosion based on thruster geometry, 
operating conditions, and magnetic field configuration. This model relies on a one-dimensional representation of the 
plasma with a fixed ionization fraction and variable ion energies based on the magnetic field distribution. Sputtering 
was modeled as the result of elastic scattering of ions by neutrals within the channel. Not all scattered ions and 
neutrals were assumed to reach the channel walls as a result of additional subsequent scattering events. 
Incorporating this phenomenon allowed the model to predict a decrease in erosion rate with time not predicted by 
only accounting for geometric effects. 

In the future this model will be used to predict the erosion of other Hall thrusters and/or operating conditions to 
determine if agreement with the SPT-100 data was simply fortuitous or if such a simplified approach can be used 
for lifetime predictions. A sensitivity analysis to the input parameters will be performed and many of the simplifying 
assumptions will be re-evaluated. As currently configured the code runs in less than two minutes suggesting that 
further increases in complexity would not necessarily be prohibitively computationally intensive. As a result we plan 
to investigate increasing the model to two-dimensions and using a plasma based model to predict neutral and plasma 
densities rather than relying on assumptions of propellant utilization fraction. 
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